{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we'll compute the solvent accessible surface area of one of the residues in our protien\n",
    "accross each frame in a MD trajectory. We're going to use our trustly alanine dipeptide trajectory for\n",
    "this calculation, but in a real-world situtation you'll probably want to use a more interesting peptide\n",
    "or protein, especially one with a Trp residue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import mdtraj as md"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use the algorithm from Shrake and Rupley for computing the SASA. Here's the function in MDTraj:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "help(md.shrake_rupley)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trajectory = md.load('ala2.h5')\n",
    "sasa = md.shrake_rupley(trajectory)\n",
    "\n",
    "print(trajectory)\n",
    "print('sasa data shape', sasa.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The computed `sasa` array contains the solvent accessible surface area for each atom in each frame of the trajectory. Let's sum over all of the atoms to get the total SASA from all of the atoms in each frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "total_sasa = sasa.sum(axis=1)\n",
    "print(total_sasa.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.pylab import *\n",
    "\n",
    "plot(trajectory.time, total_sasa)\n",
    "xlabel('Time [ps]', size=16)\n",
    "ylabel('Total SASA (nm)^2', size=16)\n",
    "show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We probably don't really have enough data do compute a meaningful [autocorrelation](http://en.wikipedia.org/wiki/Autocorrelation), but for more realistic datasets, this might be something that you want to do."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def autocorr(x):\n",
    "    \"Compute an autocorrelation with numpy\"\n",
    "    x = x - np.mean(x)\n",
    "    result = np.correlate(x, x, mode='full')\n",
    "    result = result[result.size//2:]\n",
    "    return result / result[0]\n",
    "\n",
    "semilogx(trajectory.time, autocorr(total_sasa))\n",
    "xlabel('Time [ps]', size=16)\n",
    "ylabel('SASA autocorrelation', size=16)\n",
    "show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
